Multi‐dataset identification of innovative feature genes and molecular mechanisms in keratoconus

Abstract This study aimed to identify feature genes and explore the molecular mechanisms of keratoconus (KC). We downloaded data files from NCBI GEO public database. The Limma package was used for differential expression analysis of gene profiles. Lasso regression was used to identify the feature genes. The CIBERSORT algorithm was used to infer the proportion of immune‐infiltrating cells and analyse the correlation between gene expression levels and immune cells. Related transcription factors and miRNAs of key genes were predicted using the Cistrome DB and Mircode databases. Analysis of expression differences in disease genes was based on the GeneCards database. The CMap was used to analyse targeted therapeutic drugs. IHC was performed to verify the expression levels of ATOH7 and MYRF in corneas. Exactly 593 upregulated and 473 downregulated genes were identified. Lasso regression analysis identified ATOH7, DBNDD1, RNF217‐AS1, ARL11, MYRF and SNORA74B as feature genes for KC. All key genes were correlated with immune infiltration and the levels of activated memory CD4+ T cells and plasma cells were significantly increased. miRNA, IRF and STAT families were correlated to feature genes. The expression levels of key genes were significantly correlated to KC‐related genes. Entinostat, ochratoxin‐a, diphencyprone and GSK‐3‐inhibitor‐II were predicted as potential KC medications. The expression of MYRF was significantly higher in the KC samples, contrary to the expression of ATOH7. KC is related to both immune infiltration and genetic factors. MYRF and ATOH7 were newly identified and verified feature genes of KC.

cone-shaped appearance. 3Once KC achieves to the severe stage and certifies with excessive ectasia, thinning and scarring, thereby obviously impairing visual acuity, keratoplasty is considered the last resort. 4Santodomingo-Rubido et al., illustrates the prevalence and incidence rates of KC have been approximated to be between 0.2 and 4790 per 100,000 persons and 1.5 and 25 per 100,000 persons/ year. 2 The highest prevalence and incidence rates classically find in individuals between the ages of 20 to 30 years. 5,6Thus, KC is one of the leading indications for keratoplasty worldwide and causes low vision disability. 7,8though the pathogenesis of KC remains uncertain, it is mainly considered to be related to genetic and environmental factors. 9It has been estimated that individuals with a family history of KC have a 15-67 times higher risk of developing KC than those with no family history. 10Further studies on the candidate genes of KC (SPATA7 11 and KC6 12 ) and others 13 illustrate that mutations in the presence of gene variants are required to elicit keratoconic traits. 14Traditionally, KC has been considered as a non-inflammatory corneal disease. 4wever, increasing studies have found that immune response and inflammation have important roles in the pathogenesis of KC. 15 Different kinds of immune cells and inflammatory factors have been found on the ocular surface and circulatory system of patients with KC. 16 The levels of multiple inflammatory factors in the cornea and tears of patients with KC are significantly up-regulated, and some of them (e.g.IL-10, MMP-9, TGFB1) are unique to KC. 17 Moreover, the IL-17 signalling pathway may have an important role in the pathogenesis and development of KC. 18 Recently, the rapid progress in the development of nextgeneration sequencing (NGS) technologies, including bulk RNA sequencing (bulk RNA-seq) and single-cell RNA-sequencing (scRNAseq), has shed light on complicated biological systems including cancer genomics and diverse microbial communities. 19Bulk RNA-seq is a method for the transcriptomic analysis of pooled cell populations, tissue sections or biopsies, 20 which enables comprehensive and rapid access to almost all transcript sequence information in a given tissue or organ of a species in a given state. 21scRNA-seq is used to analyse the expression profile of the cellular transcriptome at the level of a single cell, identify cell-specific markers, find rare cell types and reveal differential expression between cells.In addition, it has become an effective method to study the cellular heterogeneity of complex organisms. 22 this research, we discovered innovative feature genes of KC and explored the function of feature genes in immune infiltration and the underlying molecular mechanisms of KC by combining bulk RNA-seq and scRNA-seq.The expression levels of ATOH7 and MYRF were verified by immunohistochemical staining of corneas.
The flowchart of analysis was shown in Figure 1.

| Data download
The Series Matrix File data file of GSE204791 was downloaded from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) public database with annotation file GPL21185, which included the expression profile data of 16 patients, including 8 normal and 8 disease groups.The Series Matrix File data file of GSE77938 was downloaded from the NCBI GEO public database with the annotation file GPL18460, which included the expression profile data of 50 patients, including 25 normal and 25 disease groups.Moreover, four sample datasets of GSE146276 were downloaded from the NCBI GEO public database for singlecell analysis.
F I G U R E 1 Workflow of identification of key genes analysis in KC.DEG, different expressed genes; TF, transcription factors.

| Differential expression analysis
The Limma package is an R software package used for the differential expression analysis of gene profiles, which identifies genes with significant differential expression between different comparison groups.
We used the R package "Limma" to analyse the different molecular mechanisms of the KC data, identify differentially expressed genes between control and disease samples, select differentially expressed genes with screening criteria of p <0.05 & |log2FC| >1, and draw differential gene volcano plots and heat maps.

| Functional enrichment analysis
The metascape database (www.metas cape.org) was used to perform functional annotation of the differentially expressed genes to explore the relevance among these genes.Gene Ontology (GO) pathway and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed on differentially expressed genes.A minimum overlap of ≥3 and p ≤ 0.01 was considered statistically significant.

| Construction of the prediction model
The candidate gene set was selected, and lasso regression was applied to construct a predictive model for the diagnosis of KC.After incorporating the expression values of each specific gene, a risk score formula was constructed for each patient, and the estimated regression coefficients were weighted in the lasso regression analysis.The score for each patient was calculated according to the risk score formula, and the accuracy of the model prediction was verified by receiver operator characteristic (ROC) curve and the value of the area under ROC curve (AUC).

| Immune cell infiltration analysis
The CIBERSORT method is widely used to evaluate immune cell types and subtypes in the microenvironment.This method was according to the support vector regression principle, and deconvolution analysis was performed on the expression matrix of immune cell subtypes.It included 547 biomarkers that distinguished 22 human immune cell phenotypes, including T, B, plasma and myeloid cell subgroups.In this study, the CIBERSORT algorithm was used to analyse patient data to evaluate the relative proportions of 22 immuneinfiltrating cells and perform Spearman's correlation analysis between gene expression levels and immune cell content.

| Gene set enrichment analysis (GSEA) pathway enrichment analysis
GSEA uses predefined gene sets to rank genes according to their differential expression levels in two types of samples.It then tests whether the pre-set gene set is enriched at the top or bottom of this ranking table.In this study, GSEA was used to compare the differences in the KEGG pathway between the high-and lowexpressed samples of feature genes and explore the molecular mechanisms of key genes in the two groups of patients.The permutation number was set to 1000 and the permutation type was set to the phenotype.

| Connectivity map (CMap) drug prediction
The

| Single-cell sequencing analysis
First, the data were processed using the R package named Seurat package and the tSNE algorithm was used to analyse the positional relationships between each cluster.The clusters were annotated to cells that were important for disease development using the Celldex package.Finally, we used the logfc threshold parameter of FindAllMarkers to extract marker genes for each cell subtype from the single-cell expression profile.Genes with |avg log2FC| >1 and p < 0.05 were selected as unique marker genes for each cell subtype.For immunohistochemistry (IHC), the dissected tissues were placed in a 59°C constant temperature incubator for 60 min before dewaxing and then immersed in xylene (G20641B, GENERAL-REAGENT, NJ, US)

| Immunohistochemistry staining
and anhydrous ethanol (G73537C, GENERAL-REAGENT, NJ, US) for 10 min.The process was repeated thrice.After washing, the tissues were soaked in 200 mL of 3% H 2 O 2 and 1 mL of NaN 3 for 10 min.
Finally, the tissue sections were rinsed with anhydrous ethanol and xylene, mounted using a mounting medium (WH1161; Weiao, Shanghai, China) and covered with plastic coverslips.Images were captured using a biological microscope (CX43; Olympus, Tokyo, Japan).The positive cell area of each slide representing the protein expression levels of MYRF and ATOH7 was analysed using Image J software (https:// imagej.net/ ij/ , Fiji) under a 20x view according to the existing protocol. 24

| Statistical analysis
All statistical analyses were performed using R software (version 4.2.2).The level of statistical significance was set at p < 0.05.

| Identification of differential genes and pathway enrichment
We downloaded the GSE204791 dataset from the NCBI GEO public database, which included 16 cases: eight in the normal group and eight in the disease group.The Limma package was used to calculate the differentially expressed genes in the two groups with screening criteria of p < 0.05 & |log2FC| >1.In total, 1066 differentially expressed genes were identified, including 593 upregulated and 473 downregulated genes (Figure 2A,B), with more details provided in Material S1.Subsequently, we analysed differentially expressed genes using pathway analysis.With the Metascape database, the differentially expressed genes were primarily enriched in pathways such as secondary metabolic processes, retinoid metabolic processes and monooxygenase activity (Figure 3A).
The network of these pathways is displayed in Figure 3B.

| Identification of key genes and construction of the diagnostic model
We used the GSE204791 and GSE77938 datasets as the training and validation sets, respectively.The differentially expressed genes identified in the previous step were selected for feature selection using Lasso regression.The results showed that six genes were identified by Lasso regression as feature genes for KC, namely ATOH7, DBNDD1, RNF217-AS1, ARL11, MYRF and SNORA74B (Figure 4A-C).Except for ATOH7, the expression levels of the other key genes in the KC group were significantly upregulated (Material S1).These six genes were used as core

| Single-cell analysis of key genes
We downloaded single-cell data from GSE146276 and performed single-cell analysis using the Seurat package.The cells were clustered using the tSNE algorithm, and SingleR in the R package was used to annotate each cluster.All clusters were annotated into three cell categories: epithelial cells, keratinocytes, and fibroblasts.
Among these, epithelial cells accounted for the highest proportion (Figure 5A).The expression of key genes in these three cell types was shown in Figure 5B,C.We found that ATOH7, DBNDD1, ARL11, and MYRF were all expressed in epithelial cells and the ranking in terms of expression level was DBNDD1, ARL11, MYRF and ATOH7 (Figure 5B).ARL11 was the most expressed gene in keratinocytes and lower expression of these feature genes was observed in fibroblasts than in other cell types, except MYRF (Figure 5C).

| Immune infiltration analysis
The immunological microenvironment of the cornea, including the immune cells, cytokines, chemokines, and immune mediators, is usually altered in patients with KC. 17  immune cell infiltration and had important roles in the immune microenvironment.

| Prediction of transcription factors
We used key genes as the gene set to further explore the transcriptional factor regulatory networks involved in KC.The relevant transcription factors were predicted using the Cistrome DB online database.Among them, 99 transcription factors were predicted for DBNDD1, 81 were inferred for ARL11, 48 were predicted for MYRF and 76 were predicted for ATOH7 (Figure 9).Several immune response-related transcription factors were identified in this study.For instance, regarding the IRF family, IRF4 was found in ATOH7, DBNDD1, MYRF and ARL11 and IRF5 was predicted for MYRF (Figure 9).In the STAT family, STAT4 was found in MYRF and DBNDD1, and STAT3 was found in ARL11 (Figure 9).We constructed a comprehensive transcriptional regulatory network of key genes in KC using Cytoscape for visualization (Figure 9).Further details were provided in Material S3.

| Disease regulatory network of key genes
We obtained 1802 genes related to KC from the GeneCards database (https:// www.genec ards.org/ ).Analysis of the expression differences of disease genes revealed that the expression of genes, such as KC6, LCA5, SPATA7 and ZNF469, differed between the normal and disease groups (Figure 10A).Next, we conducted a correlation analysis between key and KC-related genes.We found that the expression levels of key genes were significantly correlated with those of KC-related genes (Figure 10B

| Upstream miRNA network analysis and potential therapeutic drug prediction
We conducted a reverse prediction of key genes using the Mircode database and obtained 49 miRNAs and 74 mRNA-miRNA relationship pairs.The results were visualized using Cytoscape (Figure 11).The data showed that miR-150 had a strong interaction with DBNDD1, SNORA74B and ARL11 (Figure 11).We divided the differentially expressed mRNAs into upregulated and downregulated groups and used the CMap database to predict drug targets for the differentially expressed genes.The results showed that the expression profiles perturbed by drugs, such as entinostat (Figure 12A), ochratoxin-a (Figure 12B), diphencyprone (Figure 12C) and GSK-3-inhibitor-II (Figure 12D), were significantly negatively correlated with disease-perturbed expression profiles, and the drugs might alleviate or even reverse the disease state.More details on miRNA and drug prediction were provided in Materials S4 & S5.

| Immunohistochemistry staining
We verified the expression levels of MYRF and ATOH7 in normal and KC corneal samples through IHC staining.The positive cell area, which represented the gene expression level, was calculated as described in the Method section.Consistent with the results mentioned above, the expression of MYRF was significantly higher in KC samples (2.13%) than in control samples (1.62%) (p < 0.01; Figure 13 A-C).In contrast, the expression of ATOH7 was significantly lower in KC samples (0.904%) than in control samples (1.22%) (p < 0.05; Figure 13D-F).Both ATOH7 and MYRF were mainly expressed in the epithelium, but more MYRF staining was observed in the stroma of the KC group (Figure 13E), in which fibroblasts accounted for a higher proportion, consistent with the results of the scRNA-seq analysis shown in Figure 5.

| DISCUSS ION
Although various therapies for KC have been developed, including wearing contact/scleral lenses, corneal cross-linking, phototherapeutic keratectomy and keratoplasty, the outcomes remain Transcriptional regulatory network of key genes in keratoconus.
unsatisfactory. 4Thus, it is essential to conduct in-depth research on the biological mechanisms and pathogenesis of KC to determine a more effective and efficient therapy for KC.
In this study, we enriched KC-related pathways using bulk RNAseq analysis of differentially expressed genes in KC.Secondary metabolic process pathways 25 are chemical reactions resulting in many of the chemical syntheses of compounds.The retinoid metabolic process 26 is associated with the synthesis of structurally similar natural RNF217-AS1, ARL11, MYRF and SNORA74B as feature genes for KC (Figure 4).Constructed from the key genes, the predictive model had an AUC of 1.0 and was verified with an AUC of 0.832, indicating that our findings provided new information regarding KC diagnosis.
The annotation results of scRNA-seq were analysed, and the key genes (ATOH7, DBNDD1, ARL11, SNORA74B and MYRF) were not differentially expressed among the cell subsets, suggesting that these genes were not specific biomarkers of the cell subsets.In terms of immune cell composition, native B cells accounted for the highest proportion with a role in the secondary adaptive response and immunological memory in human autoimmune diseases. 28Plasma cells and activated memory CD4 + T cells of the disease group were significantly higher than those of the normal group and had a positive correlation (Figure 6B,C), indicating that there was immune response in KC.CD4 + memory T cells developed by antigen exposure or the cytokine milieu communicate with various immune and non-immune cellular networks within distinct microenvironments in human chronic autoimmune diseases such as multiple sclerosis. 29milarly, plasma cells produce protective antibodies against infectious agents and pathogenic antibodies in autoimmune diseases like systemic lupus erythematosus. 30This could be explained by the abnormal expression of core genes due to the significant positive correlations between plasma cells and ARL11, DBNDD1 and MYRF (all p < 0.05; Figure 6D,F, G) and the significant positive correlations between activated memory CD4 + T cells and ARL11 and SNORA74B (both p < 0.05; Figure 6D,I) and negative correlations with ATOH7 (p < 0.01; Figure 6E).
In addition, our analysis revealed that MYRF, RNF217-AS1 and SNORA74B exhibited stronger positive correlations with immunostimulatory factors, such as CD86 (Figure 7C), in contrast to ATOH7.

These findings indicate that MYRF, RNF217-AS1 and SNORA74B
may play significant roles in promoting and maintaining immune activation in KC.CD86 is a member of the immunoglobulin superfamily that is crucial for initiating both innate and adaptive immune responses through the activation of naïve and memory T cells. 31nversely, we observed that MYRF, RNF217-AS1 and SNORA74B also showed more substantial positive correlations with immunoinhibitory factors, such as VTCN1 (Figure 7B), compared to ATOH7.
VTCN1 functions to suppress T-cell activation by inhibiting IL-2 production and influencing the cell cycles of CD4 + and CD8 + T cells. 32Interestingly, the predominant presence of plasma cells and activated memory CD4 + T cells, without significant involvement of observed in other inflammatory diseases. 33Moreover, the chemokine receptors CCR2 and CCR5, as well as F2RL1 and CXCL5, were involved in the immune regulation process of KC. 34 D'Souza, Sharon et al. also reported the significantly higher levels of inflammatory cytokine in tear fluid and implied a distinct immuno-inflammatory component in KC pathogenesis. 35The associations between KC and immune-mediated diseases have also been identified (Hashimoto's thyroiditis and inflammatory skin conditions) and confirmed (atopic conditions, including allergic rash, asthma and bronchial hyperresponsiveness and allergic rhinitis) by Claessens, Janneau L J et al. 36 Combined with the known papers and our findings, the observed changes in immune cell composition and the presence of immune responses in KC patients might originate from disruptions in the expression of these key genes.Understanding the interplay between these genes and immune modulators is crucial for elucidating the underlying mechanisms of KC and could lead to the development of targeted immunotherapeutic strategies.
To explore the biological mechanism, the GSEA results of the key genes showed an obvious trend of gene enrichment (Figure 8).
The MAPK pathway, which is enriched in MYRF, may promote inflammatory cytokines in the pathogenesis of KC. 37 The B cell receptor signalling pathway enriched in RNF217-AS1 also implied the immune response related to plasma cells existing in patients with KC. 38 Relationships between inflammation-relevant transcription factors (IRF and STAT families) and feature genes were identified in this study.IRF4 was correlated with ATOH7, DBNDD1, MYRF and ARL11 with a key regulatory function in the mature of human immune cells and the essentiality for the differentiation of T lymphocytes and B lymphocytes as well as certain myeloid cells. 39IRF5 related to MYRF acts as a molecular switch that decides whether macrophages initiate or inflammation suppress 40 (Figure 8).With regards to STAT family, STAT4, found in MYRF and DBNDD1, is essential for the development of Th1 cells from naïve CD4 + T cells 41 and STAT3, found in ARL11, is essential for the differentiation of TH17 helper T cells 42 and could be activated via the phosphorylation of serine 727 by the MAPK pathway 43 (Figure 9).Other kinds of transcription factors also related to the core genes (Figure 9), such as EBF3 44 involved in B-cell differentiation and neurogenesis; FOS 45 promoting cell migration and metastasis, and NR2F6 46 negatively regulating T cell receptors also participated in immune response.These findings implied that feature genes may alter the immune microenvironment through regulating transcription factors during KC progress.
In our analysis of the GeneCards database for genes associated with KC, we found a significant negative correlation between MYRF and SPATA7, and a significant positive correlation between MYRF and KC6 (Figure 10).This indicates that MYRF is highly expressed in the KC group, suggesting its potential role in promoting the development of KC.MYRF expression is primarily observed in mature myelinating oligodendrocytes within the central nervous system. 47Additionally, MYRF has a critical regulatory impact on immune factors and plays an important role in the immune microenvironment. 48Variants in MYRF are known to cause autosomal dominant and syndromic nanophthalmos in humans and are crucial for the development of the retina and retinal pigment epithelium (RPE). 49These findings highlight the significant role of MYRF in ocular tissue development and its potential contribution to KC.
The miRNAs of the target genes were predicted, among which miR-150 showed the strongest interaction with DBNDD1, SNORA74B and ARL11 (Figure 11).miR-150 is commonly expressed in mature B and T cells 50 and attracts and activates immune cells in breast cancer. 51The knockout of miR-150 inhibits spontaneous T cell proliferation and the progress of colitis. 52In the cornea, hsa-miR-150-5p plays a regulatory role in corneal epithelial stem cell maintenance by inhibiting the Wnt signalling pathway. 53The potential therapeutic drug predicted from key genes showed that entinostat and GSK-3-inhibitor-II had similar chemical structures (Figure 12A,C).Both Entinostat 54 and Diphencyprone 55 have been widely used in immunotherapies for many other diseases.This further illustrates that the feature genes identified in this study are potential immunotherapeutic targets for the treatment of KC.
The IHC staining of normal and KC corneas further verified our data, with significantly higher expression of MYRF in KC corneas, contrary to the expression of ATOH7 (Figure 13).We have already discussed the important role of MYRF in the immune response and the structural development of KC.ATOH7 is involved in the development of photoreceptors, retinal ganglion cells and optic nerves. 56,57tations in ATOH7 are associated with several retinal diseases, including autosomal recessive persistent hyperplastic primary vitreous 58 and familial exudative vitreoretinopathies. 59 Our findings regarding ATOH7 in KC implied an important function of ATOH7 in the normal development of the cornea.The opposite effects and expression tendencies of MYRF and ATOH7 in KC also indicated the potential role of ATOH7 in the immune regulation of KC progression.
Considering the genetic role of ATOH7 and MYRF in ocular tissues and combining with the enriched pathways in KC, the differential expression of MYRF and ATOH7 in KC may affect the synthesis of compound in KC corneas.Therefore, KC may not only be related to the immune response but also to the abnormal development or synthesis of structures owing to genetic factors.
However, due to the shortage of primary antibody types, only the protein expression levels of MYRF and ATOH7 were verified by IHC.The lack of an accurate animal model for KC also limits our further exploration in molecular mechanism of KC as well as verification of the effect of predicted drugs on KC in vivo. 60Thus, the establishment of an effective, precise and stable animal model for KC is essential.More studies investigating the regulation of immune response and specific molecular mechanisms in KC need to be accomplished in the future.
The present research has reported notable findings, but it also has limitations.Although some significant results were obtained in this study, the limited sample size still restricted the breadth and depth of this research.In addition, the identified target genes were not prospectively validated, which might affect the reliability of the results.
Future studies with more samples, larger cohorts and prospective validation are essential to be conducted in the future.Moreover, although GSK3B inhibitors have been identified as potential therapeutic targets, our study did not find the involvement of Wnt-related genes, indicating that GSK3B may operate through alternative pathways, such as those related to cellular stress responses or fibrosis. 61,62More researches into these mechanisms are also necessary.
In conclusion, we determined the feature genes of KC, constructed a diagnostic model, predicted therapeutic targets and found that the pathophysiology of KC was not only related to immune infiltration regulated by transcription factors and miRNAs correlated to feature genes but also abnormal development or synthesis of corneal structures caused by genetic factors.
Experiments on human tissues adhered to the tenets of the Declaration of Helsinki.The experimental protocols were evaluated and approved by the Ethical Committee of the Eye & ENT Hospital, Fudan University (ky2012-037).Seven KC samples collected from post-penetrating keratoplasty patients in the operation room were at stage IV (ARC (3-mm zone) <6.15, PRC (3-mm zone) <4.95, thinnest pachy <400um, best corrected distance visual acuities (BDVA) < 20/400) according to Belin ABCD Staging/Classification Parameters. 23Seven control samples were obtained from healthy corneoscleral rings from residual corneas after use at the Eye Bank of the Eye & ENT Hospital of Fudan University.
genes for subsequent experiments and a predictive model for the diagnosis of KC was constructed.The model formula was: RiskScore = ATOH7 × (−0.0638709422255986) + DBN DD1 × 0.0 069 098671 7583838 + RNF 217-AS1 × 0.00751 5606 011 362 + AR L11 × 0.00 81415 68 3882 3462 + MYRF × 0.0290 7319 99 084 615 + SN OR A74B × 0.198 60 9773078782.The results showed that the predictive model constructed from the six genes had satisfactory diagnostic performance, with an AUC of 1 (Figure4D).The GSE77938 dataset was used as the validation set to verify the predictive model externally.The results showed that the model had strong stability with an AUC of 0.832 verified by the validation sets GSE77938 (Figure4E).
Figure 7C).Furthermore, the correlations between ATOH7 and immune factors were almost contrary to MYRF, with an obvious positive correlation with immunoinhibitors (VTCN1, Pearson r = 1, p < 0.01; Figure7B), chemokines and MHC (Figure7A,B,D).These findings suggested that these key genes were closely relevant to

Next, we studied
the specific signalling pathways enriched in key genes to investigate the potential molecular mechanisms.The complete GSEA results were presented in Material S2.Highly significant pathways were selected for a centralized display.The pathways enriched in ARL11 included ascorbate and aldarate metabolism, glycine serine and threonine metabolism and glycosphingolipid biosynthesis ganglio series (Figure 8A); the pathways enriched in ATOH7 included adipocytokine signalling pathway, melanogenesis and nitrogen metabolism (Figure 8B); the pathways enriched in DBNDD1 included DNA replication, pentose and glucuronate interconversions, and chemokine signalling pathway (Figure 8C); the pathways enriched in MYRF included MAPK signalling pathway, pentose phosphate pathway and taste transduction (Figure 8D); the pathways enriched in RNF217-AS1 included B cell receptor signalling pathway, melanogenesis and olfactory transduction (Figure 8E); and the pathways enriched in SNORA74B included cadherens junction, RNA degradation and inositol phosphate metabolism (Figure 8F).

F I G U R E 3 F I G U R E 6
Pathway enrichment of differentially expressed genes based on Metascape database.(A) Primarily enriched pathways and (B) their network.F I G U R E 4 Identification of key genes and construction of predictive model.(A) Ten-fold cross validation for tuning parameter selection in LASSO model.(B) Distribution of LASSO coefficients for differentially expressed genes.(C) Bar chart of LASSO coefficients for 6 key genes.(D) ROC of predictive model for training set GSE204791 with an AUC of 1 and (E) for validation set GSE77938 with an AUC of 0.832.F I G U R E 5 Single-cell analysis of key genes.(A) Annotation of clusters into three cell categories.(B, C) Expression of key genes in cluster of epithelial cells, keratinocytes and fibroblasts.Immune infiltration analysis.(A) Proportion of immune cell content in each sample.(B) Correlations between immune cells.(C) Comparison of immune cells in normal and disease group.(D-I) Correlation between key genes and immune cells.

F I G U R E 6 Continued | 11 of 19 LYU
et al.

F I G U R E 1 0| 15 of 19 LYU
derivatives or synthetic compounds.Monooxygenase activity27 is associated with incorporation of one hydroxyl group (-OH) into substrates in many metabolic pathways and utilized by cells to metabolize arachidonic acid (i.e.eicosatetraenoic acid) into cell signalling molecules to reduce or inactivate initiated signalling molecules.It is known that the abnormal structures of epithelial and stromal tissues in KC might be relevant to the unusual synthesis of structures similar to retinol and retina.Lasso regression identified ATOH7, DBNDD1, Disease Regulatory Network of Key Genes.(A) Expression differences of disease genes between normal and disease group.(B) Correlation of expression levels between key genes and keratoconus-related genes.et al.

F I G U R E 11
Network of predicted target miRNA related to ATOH7, DBNDD1, SNORA74B and ARL11.neutrophils or macrophages, could explain the absence of classic inflammatory signs (e.g.heat, redness, swelling and pain) commonly